home *** CD-ROM | disk | FTP | other *** search
/ GFX Sensations 1 / Graphic Sensations - Volume 1.iso / tools / amiga / 3d_tools / irit40s.lha / Irit / cagd_lib / cagdcrvt.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-12-30  |  9.7 KB  |  331 lines
  1. /******************************************************************************
  2. * CagdCrvt.c - curvature computation of curves and surfaces.              *
  3. *******************************************************************************
  4. * Written by Gershon Elber, March 93.                          *
  5. ******************************************************************************/
  6.  
  7. #include "cagd_loc.h"
  8.  
  9. #define MAX_POS_REF_ITERATION 20
  10.  
  11. /******************************************************************************
  12. * Computes a scalar curve representing the curvature of a planar curve.       *
  13. * The given curve is assumed to be planar and only its x and y coordinates    *
  14. * are considered.                                  *
  15. *   Then the curvature k is equal to                          *
  16. *      .  ..    .  ..                                  *
  17. *      X  Y  -  Y  X                                  *
  18. * k =  -------------                                  *
  19. *         .2   .2  3/2                                   *
  20. *      (  X  + Y  )                                  *
  21. *   Since we cannot represent k because of the square root, we compute and    *
  22. * represent k^2.                                  *
  23. ******************************************************************************/
  24. CagdCrvStruct *CagdCrvCurvatureSqr(CagdCrvStruct *Crv)
  25. {
  26.     CagdBType
  27.         IsRational = CAGD_IS_RATIONAL_CRV(Crv);
  28.     CagdCrvStruct *Crv1W, *Crv1X, *Crv1Y, *Crv1Z,
  29.           *Crv2W, *Crv2X, *Crv2Y, *Crv2Z,
  30.           *CTmp1, *CTmp2, *CTmp3, *Numer, *Denom, *CurvatureSqr,
  31.     *Crv1Deriv = CagdCrvDerive(Crv),
  32.     *Crv2Deriv = CagdCrvDerive(Crv1Deriv);
  33.  
  34.     CagdCrvSplitScalar(Crv1Deriv, &Crv1W, &Crv1X, &Crv1Y, &Crv1Z);
  35.     CagdCrvSplitScalar(Crv2Deriv, &Crv2W, &Crv2X, &Crv2Y, &Crv2Z);
  36.     CagdCrvFree(Crv1Deriv);
  37.     CagdCrvFree(Crv2Deriv);
  38.  
  39.     CTmp1 = CagdCrvMult(Crv1X, Crv2Y);
  40.     CTmp2 = CagdCrvMult(Crv2X, Crv1Y);
  41.     CTmp3 = CagdCrvSub(CTmp1, CTmp2);
  42.     CagdCrvFree(CTmp1);
  43.     CagdCrvFree(CTmp2);
  44.     Numer = CagdCrvMult(CTmp3, CTmp3);
  45.     CagdCrvFree(CTmp3);
  46.  
  47.     CTmp1 = CagdCrvMult(Crv1X, Crv1X);
  48.     CTmp2 = CagdCrvMult(Crv1Y, Crv1Y);
  49.     CTmp3 = CagdCrvAdd(CTmp1, CTmp2);
  50.     CagdCrvFree(CTmp1);
  51.     CagdCrvFree(CTmp2);
  52.     CTmp1 = CagdCrvMult(CTmp3, CTmp3);
  53.     Denom = CagdCrvMult(CTmp1, CTmp3);
  54.     CagdCrvFree(CTmp1);
  55.     CagdCrvFree(CTmp3);
  56.     if (IsRational) {
  57.     CTmp1 = CagdCrvMult(Crv1W, Crv1W);
  58.     CTmp2 = CagdCrvMult(CTmp1, CTmp1);
  59.     CTmp3 = CagdCrvMult(CTmp2, Numer);
  60.     CagdCrvFree(CTmp1);
  61.     CagdCrvFree(CTmp2);
  62.     CagdCrvFree(Numer);
  63.     Numer = CTmp3;
  64.  
  65.     CTmp1 = CagdCrvMult(Crv2W, Crv2W);
  66.     CTmp2 = CagdCrvMult(CTmp1, Denom);
  67.     CagdCrvFree(CTmp1);
  68.     CagdCrvFree(Denom);
  69.     Denom = CTmp2;
  70.     }
  71.     if (CAGD_IS_BSPLINE_CRV(Denom)) {
  72.     CTmp1 = CagdMakePosCrvCtlPolyPos(Denom);
  73.     CagdCrvFree(Denom);
  74.     Denom = CTmp1;
  75.     }
  76.  
  77.     CagdMakeCrvsCompatible(&Denom, &Numer, TRUE, TRUE);
  78.     CurvatureSqr = CagdCrvMergeScalar(Denom, Numer, NULL, NULL);
  79.     CagdCrvFree(Denom);
  80.     CagdCrvFree(Numer);
  81.  
  82.     CagdCrvFree(Crv1X);
  83.     CagdCrvFree(Crv1Y);
  84.     CagdCrvFree(Crv2X);
  85.     CagdCrvFree(Crv2Y);
  86.     if (Crv1Z)
  87.     CagdCrvFree(Crv1Z);
  88.     if (Crv2Z)
  89.     CagdCrvFree(Crv2Z);
  90.     if (Crv1W)
  91.     CagdCrvFree(Crv1W);
  92.     if (Crv2W)
  93.     CagdCrvFree(Crv2W);
  94.  
  95.     return CurvatureSqr;
  96. }
  97.  
  98. /******************************************************************************
  99. * Computes a vectior field curve representing the curvature of a curve, in    *
  100. * the normal direction, that is kN.                          *
  101. *                .   ..      .       .   ..     .                  *
  102. *                C x C       C     ( C x C  ) x C                       *
  103. * kN = kB x N =  -----  x  ----- = --------------                  *
  104. *                  .  3      .           .  4                      *
  105. *                | C |     | C |       | C |                      *
  106. ******************************************************************************/
  107. CagdCrvStruct *CagdCrvCurvatureNormal(CagdCrvStruct *Crv)
  108. {
  109.     CagdBType
  110.         IsRational = CAGD_IS_RATIONAL_CRV(Crv);
  111.     CagdCrvStruct *CrvW, *CrvX, *CrvY, *CrvZ,
  112.     *CTmp, *Numer, *Denom, *CurvatureNormal,
  113.     *Crv1Deriv = CagdCrvDerive(Crv),
  114.     *Crv2Deriv = CagdCrvDerive(Crv1Deriv);
  115.  
  116.     CTmp = CagdCrvCrossProd(Crv1Deriv, Crv2Deriv);
  117.     CagdCrvFree(Crv2Deriv);
  118.     Numer = CagdCrvCrossProd(CTmp, Crv1Deriv);
  119.     CagdCrvFree(CTmp);
  120.     CagdCrvSplitScalar(Numer, &CrvW, &CrvX, &CrvY, &CrvZ);
  121.     CagdCrvFree(Numer);
  122.  
  123.     CTmp = CagdCrvDotProd(Crv1Deriv, Crv1Deriv);
  124.     CagdCrvFree(Crv1Deriv);
  125.     Denom = CagdCrvMult(CTmp, CTmp);
  126.  
  127.     if (IsRational) {
  128.     CagdCrvStruct *DenomCrvW, *DenomCrvX, *DenomCrvY, *DenomCrvZ;
  129.  
  130.     CagdCrvSplitScalar(Denom,
  131.                &DenomCrvW, &DenomCrvX, &DenomCrvY, &DenomCrvZ);
  132.     CagdCrvFree(Denom);
  133.  
  134.     CTmp = CagdCrvMult(CrvW, DenomCrvX);
  135.     CagdCrvFree(CrvW);
  136.     CrvW = CTmp;
  137.  
  138.     CTmp = CagdCrvMult(CrvX, DenomCrvW);
  139.     CagdCrvFree(CrvX);
  140.     CrvX = CTmp;
  141.  
  142.     CTmp = CagdCrvMult(CrvY, DenomCrvW);
  143.     CagdCrvFree(CrvY);
  144.     CrvY = CTmp;
  145.  
  146.     CTmp = CagdCrvMult(CrvZ, DenomCrvW);
  147.     CagdCrvFree(CrvZ);
  148.     CrvZ = CTmp;
  149.  
  150.     CagdCrvFree(DenomCrvW);
  151.     CagdCrvFree(DenomCrvX);
  152.  
  153.     CurvatureNormal = CagdCrvMult(CTmp, Numer);
  154.     CagdCrvFree(CTmp);
  155.     }
  156.     else {
  157.     CagdMakeCrvsCompatible(&Denom, &CrvX, TRUE, TRUE);
  158.     CagdMakeCrvsCompatible(&Denom, &CrvY, TRUE, TRUE);
  159.     CagdMakeCrvsCompatible(&Denom, &CrvZ, TRUE, TRUE);
  160.     CrvW = Denom;
  161.     }
  162.     CurvatureNormal = CagdCrvMergeScalar(CrvW, CrvX, CrvY, CrvZ);
  163.     CagdCrvFree(CrvX);
  164.     CagdCrvFree(CrvY);
  165.     CagdCrvFree(CrvZ);
  166.     CagdCrvFree(CrvW);
  167.  
  168.     return CurvatureNormal;
  169. }
  170.  
  171. /******************************************************************************
  172. * Computes a scalar curve representing the curvature sign of a planar curve.  *
  173. * The given curve is assumed to be planar and only its x and y coordinates    *
  174. * are considered.                                  *
  175. *   Then the curvature sign is equal to                          *
  176. *      .  ..    .  ..                                  *
  177. * s =  X  Y  -  Y  X                                  *
  178. ******************************************************************************/
  179. CagdCrvStruct *CagdCrvCurvatureSign(CagdCrvStruct *Crv)
  180. {
  181.     CagdBType
  182.         IsRational = CAGD_IS_RATIONAL_CRV(Crv);
  183.     CagdCrvStruct *Crv1W, *Crv1X, *Crv1Y, *Crv1Z,
  184.           *Crv2W, *Crv2X, *Crv2Y, *Crv2Z,
  185.           *CTmp1, *CTmp2, *Numer, *Denom, *CurvatureSign,
  186.     *Crv1Deriv = CagdCrvDerive(Crv),
  187.     *Crv2Deriv = CagdCrvDerive(Crv1Deriv);
  188.  
  189.     CagdCrvSplitScalar(Crv1Deriv, &Crv1W, &Crv1X, &Crv1Y, &Crv1Z);
  190.     CagdCrvSplitScalar(Crv2Deriv, &Crv2W, &Crv2X, &Crv2Y, &Crv2Z);
  191.     CagdCrvFree(Crv1Deriv);
  192.     CagdCrvFree(Crv2Deriv);
  193.  
  194.     CTmp1 = CagdCrvMult(Crv1X, Crv2Y);
  195.     CTmp2 = CagdCrvMult(Crv2X, Crv1Y);
  196.     Numer = CagdCrvSub(CTmp1, CTmp2);
  197.     CagdCrvFree(CTmp1);
  198.     CagdCrvFree(CTmp2);
  199.  
  200.     if (IsRational) {
  201.     Denom = CagdCrvMult(Crv1W, Crv2W);
  202.  
  203.     CagdMakeCrvsCompatible(&Denom, &Numer, TRUE, TRUE);
  204.     CurvatureSign = CagdCrvMergeScalar(Denom, Numer, NULL, NULL);
  205.     CagdCrvFree(Denom);
  206.     CagdCrvFree(Numer);
  207.     }
  208.     else {
  209.     CurvatureSign = Numer;
  210.     }
  211.  
  212.     CagdCrvFree(Crv1X);
  213.     CagdCrvFree(Crv1Y);
  214.     CagdCrvFree(Crv2X);
  215.     CagdCrvFree(Crv2Y);
  216.     if (Crv1Z)
  217.     CagdCrvFree(Crv1Z);
  218.     if (Crv2Z)
  219.     CagdCrvFree(Crv2Z);
  220.     if (Crv1W)
  221.     CagdCrvFree(Crv1W);
  222.     if (Crv2W)
  223.     CagdCrvFree(Crv2W);
  224.  
  225.     return CurvatureSign;
  226. }
  227.  
  228. /******************************************************************************
  229. * Given a curve that is positive, refine it until all its control points has  *
  230. * positive coefficients. Always returns a Bspline curve.                      *
  231. ******************************************************************************/
  232. CagdCrvStruct *CagdMakePosCrvCtlPolyPos(CagdCrvStruct *OrigCrv)
  233. {
  234.     int l;
  235.     CagdCrvStruct
  236.     *RefCrv = NULL;
  237.  
  238.     switch (OrigCrv -> GType) {
  239.     case CAGD_CBEZIER_TYPE:
  240.         RefCrv = CnvrtBezier2BsplineCrv(OrigCrv);
  241.         break;
  242.     case CAGD_CBSPLINE_TYPE:
  243.         RefCrv = CagdCrvCopy(OrigCrv);
  244.         break;
  245.     case CAGD_CPOWER_TYPE:
  246.     default:
  247.         FATAL_ERROR(CAGD_ERR_UNDEF_CRV);
  248.         break;
  249.     }
  250.  
  251.     for (l = 0; l < MAX_POS_REF_ITERATION; l++) {
  252.     int i, j,
  253.         Len = RefCrv -> Length,
  254.         Order = RefCrv -> Order;
  255.     CagdRType
  256.         *KV = RefCrv -> KnotVector,
  257.         *Nodes = BspKnotNodes(KV, Len + Order, Order),
  258.         *Pts = RefCrv -> Points[1];
  259.  
  260.     for (i = j = 0; i < Len; i++) {
  261.         if (Pts[i] < 0)        /* To refine at negative control points. */
  262.             Nodes[j++] = Nodes[i];
  263.     }
  264.     if (j == 0) {
  265.         IritFree(Nodes);
  266.         break;
  267.     }
  268.     else {
  269.         CagdCrvStruct
  270.             *NewCrv = CagdCrvRefineAtParams(RefCrv, FALSE, Nodes, j);
  271.  
  272.         CagdCrvFree(RefCrv);
  273.         RefCrv = NewCrv;
  274.         IritFree(Nodes);
  275.     }
  276.     }
  277.  
  278.     return RefCrv;
  279. }
  280.  
  281. /******************************************************************************
  282. * Given a planar curve, find all its inflection points by find the zero set   *
  283. * of the sign of the curvature function.                      *
  284. ******************************************************************************/
  285. CagdPtStruct *CagdCrvInflectionPts(CagdCrvStruct *Crv, CagdRType Epsilon)
  286. {
  287.     CagdCrvStruct
  288.     *CrvtrSign = CagdCrvCurvatureSign(Crv);
  289.     CagdPtStruct
  290.     *InflectionPts = CagdCrvZeroSet(CrvtrSign, 1, Epsilon);
  291.  
  292.     CagdCrvFree(CrvtrSign);
  293.  
  294.     return InflectionPts;    
  295. }
  296.  
  297. /******************************************************************************
  298. * Given a planar curve, find all its extreme curvatrue points by find the     *
  299. * extreme set of the curvature function.                      *
  300. ******************************************************************************/
  301. CagdPtStruct *CagdCrvExtremCrvtrPts(CagdCrvStruct *Crv, CagdRType Epsilon)
  302. {
  303.     if (CAGD_NUM_OF_PT_COORD(Crv -> PType) == 2) {
  304.     CagdCrvStruct
  305.         *CrvtrSqrCrv = CagdCrvCurvatureSqr(Crv);
  306.     CagdPtStruct
  307.         *ExtremCrvtrPts = CagdCrvExtremSet(CrvtrSqrCrv, 1, Epsilon);
  308.  
  309.     CagdCrvFree(CrvtrSqrCrv);
  310.  
  311.     return ExtremCrvtrPts;    
  312.     }
  313.     else if (CAGD_NUM_OF_PT_COORD(Crv -> PType) == 3) {
  314.     CagdCrvStruct
  315.         *CrvtrNormalCrv = CagdCrvCurvatureNormal(Crv),
  316.         *CrvtrMag = CagdCrvDotProd(CrvtrNormalCrv, CrvtrNormalCrv);
  317.     CagdPtStruct
  318.         *ExtremCrvtrPts = CagdCrvExtremSet(CrvtrMag, 1, Epsilon);
  319.  
  320.     CagdCrvFree(CrvtrNormalCrv);
  321.     CagdCrvFree(CrvtrMag);
  322.  
  323.     return ExtremCrvtrPts;    
  324.     }
  325.     else {
  326.     FATAL_ERROR(CAGD_ERR_ONLY_2D_OR_3D);
  327.     return NULL;
  328.     }
  329. }
  330.  
  331.